Optimal Control of Superconducting N-level quantum systems 
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Abstract. - We consider a current-biased dc SQUID in the presence of an applied time-dependent 
bias current or magnetic flux. The phase dynamics of such a Josephson device is equivalent to 
that of a quantum particle trapped in a 1— D anharmonic potential, subject to external time- 
dependent control fields, i.e. a driven multilevel quantum system. The problem of finding the 
required time-dependent control field that will steer the system from a given initial state to a 
desired final state at a specified final time is formulated in the framework of optimal control 
theory. Using the spectral filter technique, we show that the selected optimal field which induces 
a coherent population transfer between quantum states is represented by a carrier signal having a 
constant frequency but which is time-varied both in amplitude and phase. The sensitivity of the 
optimal solution to parameter perturbations is also addressed. 



Introduction. — Superconducting circuits with 
' Josephson junctions have received a lot of attention 
recently as promising candidates for scalable quantum 
bits [1,2]. An example of such a circuit is the so-called 
phase qubit [3-5], which is based on a current-biased 
, Josephson junction. The phase dynamics of this Joseph- 
' son device is analogous to that of a quantum particle 
trapped in a 1— D anharmonic potential. Preparation and 
control of the quantum states of the anharmonic well can 
; be achieved by applying time-dependent current pulses to 
the system. The device can be used as a qubit when oper- 
ated in the lowest two eigenstates of the anharmonic well. 

The energy levels beyond the lowest two can be ad- 
dressed as well. In particular, Rabi-like oscillations in the 
I multilevel limit have been observed with a current-biased 
dc-SQUID [6]. In the context of quantum information 
processing, the coupling between the computational basis 
and the states of the noncomputational subspace results in 
adverse effects on quantum gate operations [7]. However, 
there is no need to restrict to only two energy levels. A 
generalization to qudits {i.e. systems with a single particle 
Hilbert space of dimension d > 2) has been proposed for 
quantum computation [8]. In this case, the quantum infor- 
mation is encoded in higher-dimensional Hilbert spaces. 

Quantum computation requires a precise and com- 
plete control of quantum systems. The purpose of the 
present work is to apply optimal control theory to accu- 



rately transfer the populations of qudit states present in a 
current-biased dc-SQUID. The theory of optimal control 
is a well-developed field and finds numerous applications 
to the optimisation of nonlinear and highly complex dy- 
namic systems [9]. In the quantum chemistry context, 
optimal control was originally proposed by Rabitz and co- 
workers [10] as a control scheme of reaction channels and 
was extensively used in various control experiments. Op- 
timal control theory provides a systematic and flexible for- 
malism that can be used in quantum computation to gen- 
erate reliable and high precision quantum dynamics [11]. 
A very recent application deals with the optimization of a 
NOT-gate for phase qubits [12]. 

In the qudit case, population transfer can be realised 
via coherent transitions between quantum states interact- 
ing with the external control. In this letter we will use 
optimal control theory to demonstrate the possibility of a 
population transfer from the ground state to an arbitrary 
excited state of a quantum N-level system. In general, 
correlation and interference between the various pathways 
involved in the population transfer process cannot be ig- 
nored and lead to a possibly complicated dependence of 
the optimal control field on time which may be difficult 
to implement experimentally. We therefore restrict the 
frequency content of the control field using the spectral 
filter technique [13,14]. This enables us to find optimized 
control fields that may be experimentally feasible. 
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Fig. 1: (a) Current-biased DC SQUID, (b) Phase-dependent 
anharmonic potential with quantized levels. 



Model. — We consider a dc-SQUID, biased with a cur- 
rent /{,, consisting of two Josephson tunnel junctions em- 
bedded in a superconducting loop, threaded by a flux $b, 
see Fig. la. Each Josephson junction is characterized by 
its critical current Iq and capacitance Cq. Using the me- 
chanical analogy, it can be shown that the dynamics of the 
squid's phase <j> is isomorphic to that of a fictitious parti- 
cle of mass TO = 2Co (^q/Stt)^ moving in one-dimensional 
anharmonic potential (harmonic oscillator with weak cu- 
bic perturbation) [6], see Fig. lb. Here $o = h/2e denotes 
the flux quantum. Key parameters for the potential are 
the frequency of the bottom of the well ujp and the barrier 
height AU. In the presence of a time-dependent external 
magnetic flux $b(t), the quantum dynamics is described 
by the total Hamiltonian J^tot = + He where 



Ha 



—huj„ 

2 P 



P-" + X' 



He = hiUpe{t)X . 

(1) 

Here P — (^1 / y^mfkUp^ P and X = (^^Jmujp/h^ (j) are the 



reduced momentum and position operators, respectively. 
The anharmonic dimensionless coupling a can be tuned 
with the bias current; it is small compared to unity. The 
effect of a time-dependent external flux ^b{t) is included 
via the dimensionless function e{t) = ^b{t)/^o- 

The theory described by the Hamiltonian H^, the ba- 
sis Hilbert space of which is infinite, can be approximated 
by a theory described by an effective Hamiltonian (flnite 

X iV matrix) HeS = E^I;=o(V'a.|^0|V'a.)IV'm)(V'.^I the 
basis of which is flnite. Hamiltonian iJcff reproduces only 
the low-energy physics of the system [15]. With the har- 
monic oscillator eigenf unctions {ipfi) as the expansion ba- 
sis, the matrix elements of H^ can be easily calculated. 
The eigenvalues En and wave functions \n) of the effec- 
tive Hamiltonian can be obtained by diagonalizing Hcs 
using a unitary transformation Hes = O'' H^O where 



Hq — diag(ii^07 • ■ • J En-i). Once the spectrum and wave 
functions are available, all physical information can also 
be obtained, especially the matrix elements of the control 
Hamiltonian He in Eq. (1). Within the exact diagonal- 
isation of HeS, the total Hamiltonian i?tot = + He 
is transformed into the driven N-level quantum system 
Hamiltonian whose form is H = Hq + e{t)Hi with 
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dn,m\n){m\ (2) 



n,m— 



where {\n) : n = O...A^ — 1} is a complete set of 
orthonormal eigenstates, i.e. the eigenvectors of Hgg 
corresponding to the energies En. Here, the quanti- 
ties dn^m are the transition dipole moments deflned by 

dn,ni = {n\X\m) = C'nC'm /-A dxtp p_{x )xi^^{x) 

where — {n\^p) and A is on the order of lO^yh/mujp. 
The off-diagonal elements of Hi induce transitions between 
energy eigenstates; the diagonal elements renormalize the 
energy eigenvalues, an effect known as the Stark shift. 

Quantum optimal control problem. Let time 
t be in the interval [tj = 0,ti?], for time tp fixed. An 
arbitrary state of the system at time t can be represented 
by the density matrix p acting on , the Hilbert space 
of dimension N. The density matrix evolves according to 
the Liouville-von Neumann equation 



ihp = H,p ; p{0) = pi 



(3) 



where pj is the initial state of the system. We will 
use optimal control theory to design a control field e{t) 
which drives our system from an initial state pi at time 
i/ = to a desired target state at specified final time 
tp. The problem can be formulated in terms of a cost 
functional that also takes into account experimental con- 
straints. Minimizing this cost functional leads to the de- 
sired physical target, thereby satisfying the constraints. 

The question as to whether or not there exists a control 
that steers the system to a given goal is of crucial impor- 
tance. For iV-level quantum systems subject to a single 
control, this question has been addressed in [16,17]. In 
the absence of dissipative effects, the answer is affirmative 
because the dynamical Lie group generated by iHo and 
iHi is isomorphic to the unitary group U{N), which is 
compact. A target state pp can be dynamically reached 
from Pi if there exists a unitary operator U S U{N) such 
that Pp — ij{tp)piU{tp)] [17]. Actually U{tp) represents 
the time evolution operator obeying itself the Schrodinger 
equation with the initial condition C/(0) = 1. The for- 
mal solution to the Schrodinger equation can be written 



as[/(iF)=r|exp[-i/J^ 
denotes time ordering. 



drH 



where the symbol T 



Pontryagin minimum principle. Suppose the 
system is prepared at time t/ = in the initial state 



p-2 



pi. The objective is to compute an appropriate time- 

dcpcndont control function e(f) steering the system from 
the initial state pi into a target state pp sX fixed final time 
tp- The corresponding cost functional may be written as 



dt (4) 



where is the Frobenius norm: H^Hl' = Tr^'I'A = 

Sij l^y l^- Here, the first term represents the deviation 
between the state of the system at final time p{tF) and 
the target state pp, whereas the second integral term pe- 
nalizes the field fluency with a generally time-dependent 
weight a. We will illustrate the physical meaning of a{t) 
for a specific example below. Minimizing the first term 
{i.e. the error) is equivalent to maximizing the state trans- 
fer fidelity F = Tr {p{tp)pp}. Our overall task is to find 
the control e{t) that minimizes J[e{t)] and satisfies both 
the dynamic constraint and the boundary condition (3). 
An optimal solution of this problem can be obtained us- 
ing the first order optimality conditions in the form of the 
Pontryagin minimum principle (PMP) [18,19]. These con- 
ditions arc formulated using a scalar pseudo Hamiltonian 
which may in the present case be cast in the form 



n{p,s,X) := -as'^ + Tr 



A 

ih 



H,p 



(5) 



where the adjoint state variable A is an operator Lagrange 
multiplier introduced to implement the constraint (3). 
The PMP states that the necessary conditions of optimal- 
ity for the control problem defined above are as follows: 



A 



in 



H,p 



m = pi; (6) 

X{tp) = p{tp) - pp; (7) 



= dsH = a£ -I- Im Tr 



A 



Hi,p 



(8) 



The last condition is also equivalent to a vanishing gradi- 
ent of the functional J with respect to the control e, 



5J 



= Re r den{p{t),\{t),e{t)) 
Jo 



dt 



(a{t)s{t)+lmTi ^ ^ 



Hi,p{t) 



(9) 

dt. (10) 



Numerically, an iterative procedure based on successive 

linearization must be employed to find the optimal control. 
Here we will use a gradient-based method in order to find 
a solution to the system of the necessary conditions of 
optimality. Eqs. (6)- (8). More precisely, we have used the 
L-BFGS-B routine which is based on a bound constraint 
quasi- Newton method with BFGS update rule [20]. This 
routine is appropriate and efficient for solving constrained 
as well as unconstrained problems. 



Population transfer. — We now drive the system 

from the groimd state [0) into one of the excited states 
I n), for n = 1, 2, 3, ... as an illustrative example of the effi- 
ciency of the control field generated by the optimal control 
algorithm. Suppose the system is in state pj = |0)(0| at 
time t = 0. The objective is to force the system to state 
Pp = \n){n\ for given n at time tp. 

A variety of experimental constraints may be imposed 
in an optimal control problem in order to select control 
fields that are feasible from a practical point of view [14]. 
The purpose of the multiplier a(t) in the cost func:tional 
defined in Eq. (4) is to force the control field to approach 
zero at the initial and final time in accordance with the 
experiment. For this we use the shape function [18] 

ait) =ao + ai (exp [-t/r] + exp [-{tp - t)/T]) , (11) 

where the positive constants aj are the penalty parameters 
and T is a rise time. The role of ao is to penalize high 
control field values throughout the time interval [0,ii?]; 
«! together with the exponential terms enforces the field 
to be nearly zero at the boundaries of the interval [0,tF] 
while simultaneously turning on and off the field smoothly. 

Because all the matrix elements of the interaction 
Hamiltonian Hi defined in Eq. (2) are different from zero 
the population transfer, for example from the ground state 
|0) to the excited state |4), will involve several path- 
ways. The population transfer can be realized via a di- 
rect transition: |0) — > |4) or via indirect transitions: 
|0) |1) |2) |3) |4), |0) — |6) ^ |4) . . . 
As a result, the time structure of the control field emerg- 
ing from optimal control theory will generally be compli- 
cated, due to quantum mechanical interference between 
the pathways the field employs in the process |0) — > |4). 
The resulting control field is generally characterized by a 
frequency content that is not readily implemented exper- 
imentally. In order to reduce the control field complexity, 
we resort to the spectral filter technique. A convenient 
way to restrict the optimal field to a single desired fre- 
quency u)Q is to filter the gradient by the formula [13, 14], 
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where 
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5e{t) 
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\/27r J-oo 
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fc(i)' 



'dt 



(12) 



(13) 



is the Fourier transform of the gradient of the cost func- 
tional defined in Eq. (10) and 



g{u>) = e-T^'^-'^o)' + e-T('^+'^o)' 



(14) 



is a Gaussian frequency filter centered around ±ojo. Here 
7 is positive constant, large compared to the unity so that 

g{uj) = for Lu ^ ±a;o and g{Lu) = 1 for lu = ztcjo- During 
the optimization process, the control variable is updated 
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by a filtered gradient in each iteration. Thus, every spec- 
tral component in the control variable is eliminated except 
the components around uj = iwo- The resulting optimized 
control field has a simple intuitive interpretation. 

The restriction of the spectrum to a single frequency 
LOo = (El — Eo)/h simplifies the time structure of the opti- 
mal control. Specifically, the time dependence of the con- 
trol field can be interpreted as an input signal represented 
by a single oscillation of the form 

Soptit) = A{t) cos [ujot + ip{t)] (15) 

where the amplitude A{t) and phase ip{t) vary slowly with 
time compared to luq. The goal of complex demodula- 
tion is to extract the amplitude and phase as function 
of time [21]. First, the original signal eopt(i) is multi- 
plied by a complex modulation of frequency ujq, yielding 
f{t) = eoptl^e""^"* = ^(i)[e'^(*) _^ g-2iu.ot-«v(t)]/2. Pass- 
ing the resulting signal f{t) trough an ideal low-pass filter 
of cutoff frequency lUc < t^o leads to g{t) = A(t)e"^'^*V2 = 
a(t) + if3{t). The time-dependent amplitude and phase 
can then be calculated as A{t) = 2y^a{t)^ + P{t)^ and 
ip{t) = arctan [/3(t)/a(i)], respectively. 

Results. — We now illustrate the above concepts with 
the aid of some representative examples. Throughout 
this section, for the numerical simulations of population 
transfer from the ground state |0) to one of the excited 
states |n), we use the following nominal parameters, typi- 
cal for experiments with current-biased SQUIDs [6] : an- 
harmonicity a = 0.0325; target time tp = 500/ijjp, chosen 
to be short enough to avoid substantial relaxation and 
decoherence phenomena; time step At — 3.0 x 10^^/wp 
corresponding to M — 2^^ as the number of mesh points; 
frequency filter luq — {Ei — Eo)/h; penalties factors ao — 
10~^; ai = 10^ and rise time t = 10^/wp. 

We start by considering the population transfer in a 
two- level system, N — 2. It is well-known that popu- 
lation transfer can be obtained within the rotating wave 
approximation (RWA) by applying a so-called 7r-pulse: a 
time-dependent signal of frequency loo = Ei — Eq (i.e. res- 
onant with the level spacing) with duration equal to half 
of the so-called Rabi period. It is useful to study this sim- 
ple example with the help of optimal control theory, as it 
enables one to compare to the known result for a 7r-pulse 
and hence test the numerical implementation; it also sheds 
some light on the functioning of the control procedure. 

We first study the response of a two-level system to a 
TT-pulse. The pulse is shown Fig. 2b, its complex demod- 
ulation is shown in panels (c) and (d). The response of 
the two-level system is obtained numerically and shown in 
Fig. 2a; we see that perfect population transfer is achieved. 
We subsequently use the 7r-pulse of Fig. 2b as a guess for 
an optimal control simulation, the results of which are 
shown in Fig. 3. It can be seen that the application of op- 
timal control slightly modifies the original 7r-pulse. Specif- 
ically, the overall amplitude is smaller as a result of the 
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Fig. 2: Population transfer from the ground state |0) to the 
excited state following a 7r-pulse. Panel (a) shows the 
numerically obtained evolution of populations, (b) and (c), re- 
spectively, show the TT-pulse and its numerically obtained power 
spectrum. The amplitude of the pulse versus time is also shown 
in panel (b) while the numerically obtained time evolution of 
the phase is displayed in the panel (d). 
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Fig. 3: Optimal control of population transfer from the ground 
state |0) to the excited state (a) shows the evolution of 
populations, (b) and (c), respectively, show the selected control 
field and its power spectrum. The amplitude of the control 
field versus time is also shown in the panel (b) while the time 
evolution of the phase is displayed in the panel (d). 

minimal amplitude constraint we imposed; the pulse du- 
ration has increased in order to conserve the total area of 
the envelope and hence the tt nature of the pulse. 

We next wish to implement a population transfer from 
the ground state |0) to the excited state |1) in a multi- 
level system with = 6. As an initial guess we use the 
TT-pulse of Fig. 2b. The response of the seven-level sys- 
tem to this guess is shown in Fig. 4a: the presence of the 
levels beyond n = 1 clearly leads to strong contamination 
effects that limit the efficiency of the transfer. We then 
use the initial guess as a starting point for an optimal 
control simulation; the resulting optimal pulse is shown in 
Fig. 5b, the result of the corresponding population trans- 
fer in panel (a). In order to avoid contamination, both the 
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Fig. 4: Population transfer from the ground state |0) to the Fig. 6: Control of opulation transfer from the ground state |0) 



excited state |1) in a seven-level system, using a vr-pulse; (a) 
(b), (c) and (d), as in the previous Figure. 



to the excited state |4) in a seven-level system; (a), (b), (c) and 
(d), as in the previous Figure. 
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Fig. 5: Control of population transfer from the ground state 
|0) to the excited state |1) in a seven-level system, using the vr- 
pulse as a gues; (a), (b), (c) and (d), as in the previous Figure. 



pulse amplitude and its phase vary smoothly as a function 
of time, as shown in panels (c) and (d). We verified that 
a transfer from |0) to |1) can also be achieved in principle 
with a TT-pulse similar to the one shown in Fig. 4b. How- 
ever, the overall amplitude of the required pulse should 
be much smaller (by about a factor ten) in order to avoid 
non-resonant transitions to the higher levels. Accordingly 
the resulting duration would be ten times longer, thereby 
exceeding typical relaxation and decoherence times. 

Finally, in Fig. 6, we show an optimal control simulation 
for the population transfer from the ground state |0) to the 
excited state |4) . As can be seen by comparing panel (a) of 
Fig. 6 with that of Fig. 5, this transfer involves enhanced 
occupation during manipulation of the excited states of 
the system. We therefore expect transfers to higher states 
to be less robust against noise than transfers to low-lying 
excited states (see also the next Section). 
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Fig. 7: Sensitivity of the population transfer to a small ampli- 
tude perturbation of the optimal control field. 



Sensitivity analysis. — In applications, the system 
parameters are usually not fixed but may be subject to 
perturbations and noises. For instance, the environment 
of the dc-SDUID induces time-dependent fluctuations of 
the bias current and flux. Because of the external per- 
turbations, practical devices are not capable of operating 
precisely neither at the prescribed system parameters nor 
at the computed control field. Then, it is of great impor- 
tance to know the sensitivity of the optimal solution with 
respect to perturbations of any of the system parameters. 
In this section, we shall only consider the influence of slow, 
adiabatic fluctuations. In the presence of this so-called 
adiabatic noise, the system parameters remain constant 
during a given manipulation but fluctuate during repeti- 
tive measurements needed to obtain quantum statistics. 

In order to develop some feeling for the sensitivity of 
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Fig. 8: Sensitivity of the population transfer fidelity to the 
random fiuctuations of the anharmonic couphng. 

the model system discussed here with respect to adiabatic 
noise, we show in Fig. 7 the effect of a small static am- 
plitude fluctuation Ae (of about 1%) with respect to the 
optimal amplitude on the result of the optimal control of a 
seven- level system when a population transfer from |0) to 
|1) is sought [panel (a) and (b)] or when a transfer from |0) 
to |4) is sought [panel (c) and (d)]. As is to be expected, 
the sensitivity to such a fluctuation is larger in the latter 
case, which involves more excited states during transfer. 

We finally consider the effect of a slow variation of a 
system parameter during application of an optimal con- 
trol pulse. Specifically, we assume variations of the anhar- 
monic coupling a to occur when the same optimal con- 
trol sequence is repeated many times. Because the anhar- 
monic coupling depends on the so-called working point, 
every random perturbation that changes the bias current 
or flux changes the parameter a leading to uncertainties 
in the eigenvalues En and the transition dipole moments 
dn,m iu Eq. (2). In Fig. 8 we show a plot of the state 
transfer fidelity F as a function of cr, calculated for a con- 
trol field optimized for the case a — W — 0.0325. The 
fidelity drops as a deviates from a; we see again that it 
decreases faster as transitions to higher levels are consid- 
ered. From this result we can calculate the average fidelity 
F = J^Y.iLlF^ = igEf=iTr{p»(tF)pd}, when the op- 
timized control sequence is repeated M times, assuming 
(7 to be a normally distributed random variable with a 
mean ct — 0.0325 and a standard deviation Act = ct/16. 
This choice corresponds to a typical experimental low- 
frequency noise present in phase qubits [4,5], yielding a 
Q-factor of about 1000. We find the following average 
state transfer fidelities: F|o)^|i) = 85%, F|o)^|2) = 73%, 
F|o)^|3) = 60%, and F|o)^|4) = 55%. If the low-frequency 
noise will be reduced only by a factor five (Q-factors 
^ 5000), all the above fidelities remain larger than 95%. 

Conclusion. — Using optimal control theory, we have 
studied the possibility of a population transfer from the 
ground state to an arbitrary excited state of a supercon- 
ducting quantum N-level system. We have found that 
such state transfer can be obtained with good fidelity, us- 



ing optimized pulses which can be realized using existing 
microwave technology. We have considered the effects of 
low-frequency noise and found that it reduces the average 
state transfer fidelity. Using parameters describing actual 
phase qubits, we find a loss of fidelity of about 15 % for 
the transfer from level |0) to |1) up to 45 %, for a transfer 
from level |0) to |4). However, a substantial improvement 
of fidelity is possible by slightly decreasing the effects of 
low-frequency noise. 

We thank R. Fazio and F. Taddei for useful discus- 
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